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Abstract. La Buddc's method computes the characteristic polynomial of a real matrix A in 
two stages: first it applies orthogonal similarity transformations to reduce A to upper Hessenberg 
form H, and second it computes the characteristic polynomial of H from characteristic polynomials 
of leading principal submatrices of H . If A is symmetric, then H is symmetric tridiagonal, and La 
Budde's method simplifies to the Sturm sequence method. If A is diagonal then La Budde's method 
reduces to the Summation Algorithm, a Horner-like scheme used by the MATLAB function poly to 
compute characteristic polynomials from eigenvalues. 

We present recursions to compute the individual coefficients of the characteristic polynomial 
in the second stage of La Budde's method, and derive running error bounds for symmetric and 
nonsymmetric matrices. We also show that La Budde's method can be more accurate than poly, 
especially for indefinite and nonsymmetric matrices A. Unlike poly, La Budde's method is not 
affected by illconditioning of eigenvalues, requires only real arithmetic, and allows the computation 
of individual coefficients. 
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1. Introduction. We present a little known numerical method for computing 
characteristic polynomials of real matrices. The characteristic polynomial of a n x n 
real matrix A is defined as 

p(A) = det(A7 - A) = A™ + ciA"" 1 + • • • + c„_iA + c„, 

where I is the identity matrix, c\ — — trace(A) and c„ = (— f )" det(A). 

The method was first introduced in 1956 by Wallace Givens at the third High 
Speed Computer Conference at Louisiana State University [9... According to Givens, 
the method was brought to his attention by his coder Donald La Budde p 302]. 
Finding no earlier reference to this method, we credit its development to La Budde 
and thus name it "La Budde's method" . 

La Budde's method consists of two stages: In the first stage it reduces A to upper 
Hessenberg form H with orthogonal similarity transformations, and in the second 
stage it computes the characteristic polynomial of H . The latter is done by computing 
characteristic polynomials of leading principal submatrices of successively larger order. 
Because H and A are similar, they have the same characteristic polynomials. If A 
is symmetric, then H is a symmetric tridiagonal matrix, and La Budde's method 
simplifies to the Sturm sequence method [8]. If A is diagonal then La Budde's method 
reduces to the Summation Algorithm, a Horner-like scheme that is used to compute 
characteristic polynomials from eigenvalues [27] . The Summation Algorithm is the 
basis for MATLAB's poly command, which computes the characteristic polynomial 
by applying the Summation Algorithm to eigenvalues computed with eig. 

We present recursions to compute the individual coefficients of the characteristic 
polynomial in the second stage of La Budde's method. La Budde's method has a 
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number of advantages over poly. First, a Householder reduction of A to Hessenberg 
form H in the first stage is numerically stable, and it does not change the condition 
numbers [19 of the coefficients c k with respect to changes in the matrix. In contrast 
to poly, La Budde's method is not affected by potential illconditioning of eigenvalues. 
Second, La Budde's method allows the computation of individual coefficients c k (in 
the process, c±, . . . ,c k - X are computed as well) and is substantially faster if k <C n. 
This is important in the context of our quantum physics application, where only a 
small number of coefficients are required (TOJ §1], [2"Tl [2"2"] . 

Third, La Budde's method is efficient, requiring only about 5n 3 floating point 
operations and real arithmetic. This is in contrast to poly which requires complex 
arithmetic when a real matrix has complex eigenvalues. Most importantly, La Budde's 
method can often be more accurate than poly, and can even compute coefficients of 
symmetric matrices to high relative accuracy. Unfortunately we have not been able 
to derive error bounds that are tight enough to predict this accuracy. 

In this paper we assume that the matrices are real. Error bounds for complex 
matrices are derived in [26, §6]. 

Overview. After reviewing existing numerical methods for computing charac- 
teristic polynomials in we introduce La Budde's method in $3] Then we present 
recursions for the second stage of La Budde's method and running error bounds, for 
symmetric matrices in 2] and for nonsymmetric matrices in In ^BJ we present 
running error bounds for both stages of La Budde's method. We end with numerical 
experiments in <j7]that compare La Budde's method to MATLAB's poly function and 
demonstrate the accuracy of La Budde's method. 

2. Existing Numerical Methods. In the nineteenth century and in the first 
half of the twentieth century characteristic polynomials were often computed as a 
precursor to an eigenvalue computation. In the second half of the twentieth century, 
however, Wilkinson and others demonstrated that computing eigenvalues as roots 
of characteristic polynomials is numerically unstable |31[ 132] . As a consequence, 
characteristic polynomials and methods for computing them fell out of favor with the 
numerical linear algebra community. We give a brief overview of these methods. They 
can be found in the books by Faddeeva [5J, Gantmacher [5J, and Householder [T7] , 

2.1. Leverrier's Method. The first practical method for computing charac- 
teristic polynomials was developed by Leverrier in 1840. It is based on Newton's 
identities pSJ (7.19.2)] 

ci = -trace(A), c k = -- trace (A k + c x A k ~ x -\ h c k _ x A) , 2 < k < n. 

k 

The Newton identities can be expressed recursively as 

c k = -- tracc(AB fc _i), where B x = A + c x I, B k = AB k _ x + c k I. 

Leverrier's method and modifications of it have been rediscovered by Faddeev and 
Sominskh, Frame, Souriau, and Wegner, see [18], and also Horst [15]. Although Lev- 
errier's method is expensive, with an operation count proportional to n 4 , it continues 
to attract attention. It has been proposed for computing A" 1 , sequentially [3] and 
in parallel [4]. Recent papers have focused on different derivations of the method 
[TJ [Tni US] , combinatorial aspects [53] , properties of the adjoint [T3] , and expressions 
of p(X) in specific polynomial bases [2], 
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In addition to its large operation count, Leverrier's method is also numerically 
unstable. Wilkinson remarks [3TJ §7.19]: 

"We find that it is common for severe cancellation to take place when 
the Ci are computed, as can be verified by estimating the orders of 
magnitudes of the various contributions to Cj." 
Wilkinson identified two factors that are responsible for the numerical instability of 
computing c^: errors in the computation of the trace, and errors in the previously 
computed coefficients c\, . . . , Ck-i- 

Our numerical experiments on many test matrices corroborate Wilkinson's obser- 
vations. We found that Leverrier's method gives inaccurate results even for coefficients 
that are well conditioned. For instance, consider the n x n matrix A of all ones. Its 
characteristic polynomial is p(X) — A™ — nA n_1 , so that C2 = • • • = c„ = 0. Since A 
has only a single nonzero singular value u\ = n, the coefficients Ck are well conditioned 
(because n — 1 singular values are zero, the first order condition numbers with respect 
to absolute changes in the matrix are zero [El Corollary 3.9]). However for n = 40, 
Leverrier's method computes values for C22 through C40 in the range of 10 18 to 10 47 . 

The remaining methods described below have operation counts proportional to 

n 3 . 

2.2. Krylov's Method and Variants. In 1931 Krylov presented a method 
that implicitly tries to reduce A to a companion matrix, whose last column contains 
the coefficients of p(X). Explicitly, the method constructs a matrix K from what 
we now call Krylov vectors: v, Av, A 2 v, . . . where v ^ is an arbitrary vector. Let 
m > 1 be the grade of the vector, that is the smallest index for which the vectors 
v, Av, . . . A m ~ 1 v are linearly independent, but the inclusion of one more vector A m v 
makes the vectors linearly dependent. Then the linear system 

Kx + A m v = 0, where K = (v Av ... A m ~ 1 v) 

has the unique solution x. Krylov's method solves this linear system Kx — —A m v for 
x. In the fortunate case when m = n the solution x contains the coefficients of p(A), 
and Xi = —c n ^i + \. If to < n then x contains only coefficients of a divisor of p(X). 

The methods by Danilevskh, Weber- Voetter, and Bryan can be viewed as partic- 
ular implementations of Krylov's method |171 §6], as can the method by Samuelson 
[28]. 

Although Krylov's method is quite general, it has a number of shortcomings. 
First Krylov vectors tend to become linearly dependent, so that the linear system 
Kx = —A m v tends to be highly illconditioned. Second, we do not know in advance 
the grade to of the initial vector v; therefore, we may end up only with a divisor of 
p(A). If A is derogatory, i.e. some eigenvalues of A have geometric multiplicity 2 or 
larger, then every starting vector v has grade m < n, and Krylov's method does not 
produce the characteristic polynomial of A. If A is non derogatory, then it is similar to 
its companion matrix, and almost every starting vector should give the characteristic 
polynomial. Still it is possible to start with a vector v of grade m < n, where Krylov's 
method fails to produce p(X) for a non derogatory matrix A [Til Example 4.2]. 

The problem with Krylov's method, as well as the related methods by Danilevskii, 
Weber- Voetter, Samuelson, Ryan and Horst is that they try to compute, either im- 
plicitly or explicitly, a similarity transformation to a companion matrix. However, 
such a transformation only exists if A is nonderogatory, and it can be numerically 
stable only if A is far from derogatory. It is therefore not clear that remedies like 
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those proposed for Danilevskii's method in [TJ], [TBI p 36], (55], [SB §7-55] would be 
fruitful. 

The analogue of the companion form for derogatory matrices is the Frobenius 
normal form. This is a similarity transformation to block triangular form, where 
the diagonal blocks are companion matrices. Computing Frobenius normal forms is 
common in computer algebra and symbolic computations e.g. [7J, but is numerically 
not viable because it requires information about Jordan structure and is thus an 
illposed problem. This is true also of Wiedemann's algorithm [501 130] ; which works 
with u T A z b, where u is a vector, and can be considered a "scalar version" of Krylov's 
method. 

2.3. Hyman's method. Hyman's method computes the characteristic polyno- 
mial for Hessenberg matrices [3TJ §7.11]. The basic idea can be described as follows. 
Let B be a n x n matrix, and partition 



D 



n- 1 1 
I ( b\ b 12 



n-1 ft b 



If B 2 is nonsingular then det(B) = (-l)™" 1 det(S 2 )(6i 2 - b{B^ 1 b 2 ). Specifically, if 
B — XI — H where H is an unreduced upper Hessenberg matrix then B 2 is nonsingular 
upper triangular, so that det(£?2) = (— l) n ~ 1 h 2 i • ■ • h n , n -i is just the product of the 
subdiagonal elements. Thus 

p(X) =h 2 i--- hn >n -i(bi2 - b[B^ x b 2 ). 

The quantity B^~ x b 2 can be computed as the solution of a triangular system. However 
b\, B 2 , and b 2 are functions of A. To recover the coefficients of X 1 requires the solution 
of n upper triangular systems [24] . 

A structured backward error bound under certain conditions has been derived in 
[24) . and iterative refinement is suggested for improving backward accuracy. However, 
it is not clear that this will help in general. The numerical stability of Hyman's method 
depends on the condition number with respect to inversion of the triangular matrix 
B 2 . Since the diagonal elements of B 2 are h 2 \, . . . , /i„. n _i, B 2 can be ill conditioned 
with respect to inversion if H has small subdiagonal elements. 

2.4. Computing Characteristic Polynomials from Eigenvalues. An ob- 
vious way to compute the coefficients of the characteristic polynomial is to compute 
the eigenvalues Xi and then multiply out 11™= l — ^»)' The MATLAB function 
poly does this. It first computes the eigenvalues with eig and then uses a Horner-likc 
scheme, the so-called Summation Algorithm, to determine the Ck from the eigenvalues 
Xj as follows: 

c = [1 zeros(l,n)] 

for j = l:n 

c(2:(j+D) = c(2:(j + D) - Ay.*c(l:j) 
end 

The accuracy of poly is highly dependent on the accuracy with which the eigen- 
values Xj are computed. In [27 , §2.3] we present perturbation bounds for characteristic 
polynomials with regard to changes in the eigenvalues, and show that the errors in 
the eigenvalues are amplified by elementary symmetric functions in the absolute val- 
ues of the eigenvalues. Since eigenvalues of non-normal (or nonsymmetric) matrices 
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are much more sensitive than eigenvalues of normal matrices and are computed to 
much lower accuracy, poly in turn tends to compute characteristic polynomials of 
non-normal matrices to much lower accuracy. As a consequence, poly gives useful 
results only for the limited class of matrices with wellconditioned eigenvalues. 

3. La Budde's Method. La Budde's method works in two stages. In the first 
stage it reduces a real matrix A to upper Hessenberg form H by orthogonal similarity 
transformations. In the second stage it determines the characteristic polynomial of H 
by successively computing characteristic polynomials of leading principal submatrices 
of H . Because H and A are similar, they have the same characteristic polynomials. 
If A is symmetric, then H is a symmetric tridiagonal matrix, and La Budde's method 
simplifies to the Sturm sequence method. The Sturm sequence method was used by 
Givens [8] to compute eigenvalues of a symmetric tridiagonal matrix T, and is the 



basis for the bisection method [10l §§8.5.1, 8.5.2]. 

Givens said about La Budde's method [HI p 302]: 

Since no division occurs in this second stage of the computation and 
the detailed examination of the first stage for the symmetric case [...] 
was successful in guaranteeing its accuracy there, one may hope that 
the proposed method of getting the characteristic equation will often 
yield accurate results. It is, however, probable that cancellations of 
large numbers will sometimes occur in the floating point additions 
and will thus lead to excessive errors. 
Wilkinson also preferred La Budde's method to computing the Frobenius form. 
He states [3U §6.57]: 

We have described the determination of the Frobenius form in terms 
of similarity transformations for the sake of consistency and in order 
to demonstrate its relation to Danilewski's method. However, since 
we will usually use higher precision arithmetic in the reduction to 
Frobenius form than in the reduction to Hessenberg form, the reduced 
matrices arising in the derivation of the former cannot be overwritten 
in the registers occupied by the Hessenberg matrix. 
It is more straightforward to think in terms of a direct derivation of 
the characteristic polynomial of H . This polynomial may be obtained 
by recurrence relations in which we determine successively the char- 
acteristic polynomials of each of the leading principal submatrices H r 



No special difficulties arise if some of the [subdiagonal entries of H] 
are small or even zero. 

La Budde's method has several attractive features. First, a Householder reduction 
of A to Hessenberg form H in the first stage is numerically stable jTUJ §7.4.3], [3~Tl 
§6.6]. Since orthogonal transformations do not change the singular values, and the 
condition numbers of the coefficients Ck to changes in the matrix are functions of 
singular values [19] , the sensitivity of the Ck does not change in the reduction from A 
to H . In contrast to the eigenvalue based method in §2.41 La Budde's method is not 
affected by the conditioning of the eigenvalues. 

Second, La Budde's method allows the computation of individual coefficients Ck 
(in the process, c\, . . . , Ck-i are computed as well) and is substantially faster if k <C n. 
This is important in the context of our quantum physics application, where only a 
small number of coefficients are required [El §1], [211 122j. 

Third, La Budde's method is efficient. The Householder reduction to Hessenberg 
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(r = l,...,n) of H. [... 



form requires 10n 3 /3 floating point operations [TO] §7.4.3], while the second stage 
requires n 3 /6 floating point operations [3TJ §6.57] (or 4n 3 /3 flops if A is symmetric 
[T01 §8.3.1]). If the matrix A is real, then only real arithmetic is needed - in contrast 
to eigenvalue based methods which require complex arithmetic if a real matrix has 
complex eigenvalues. 

4. Symmetric Matrices. In the first stage, La Budde's method reduces a real 
symmetric matrix A to tridiagonal form T by orthogonal similarity transformations. 
The second stage, where it computes the coefficients of the characteristic polynomial 
of T, amounts to the Sturm sequence method [8j. We present recursions to compute 
individual coefficients in the second stage of La Budde's method in N4.ll describe our 
assumptions for the floating point analysis in HA.21 and derive running error bounds 
in ^3"! 

4.1. The Algorithm. We present an implementation of the second stage of La 
Budde's method for symmetric matrices. Let 

fax ft \ 

ft ol 2 ft 

T= '■• '■• '■• 

\ ft. OL n ) 

be a n x n real symmetric tridiagonal matrix with characteristic polynomial p(A) = 
dct(A7 — T). In the process of computing p(X), the Sturm sequence method computes 
characteristic polynomials Pi{\) = det(A7 — T{) of all leading principal submatrices Ti 
of order i, where p n {X) = p(X). The recursion for computing p(X) is [8], [10] (8.5.2)] 

p (A) = 1, pi(X) = A - at 

Pi (X) = (X - ai)pi-i(X) -/3? Pi _ 2 (A), 2<i<n. (4.1) 

In order to recover individual coefficients of p(X) from the recursion (|4.1[) . we identify 
the polynomial coefficients 



and 



p(X) = A" + ciA™ 1 H h Cn-iX + c„ 



Pt(X) = A' + c^X 1 - 1 + ■■■ + cj^A + cf , 1 < i < n, 



where ejj. = Ck- Equating like powers of A on both sides of (|4.1j) gives recursions 
for individual coefficients Ck, which are presented as Algorithm 1. In the process, 
ci, . . . , Ck-i are also computed. 

If T is a diagonal matrix then Algorithm 1 reduces to the Summation Algorithm 
[27l Algorithm 1] for computing characteristic polynomials from eigenvalues. The 
Summation Algorithm is the basis for MATLAB's poly function, which applies it 
to eigenvalues computed by eig. The example in Figure 14.11 shows the coefficients 
computed by Algorithm 1 when n = 5 and k = 3. 



6 



Algorithm 1 La Budde's method for symmetric tridiagonal matrices 
Input: n x n real symmetric tridiagonal matrix T, index k 
Output: Coefficient ci, . . . , c* of p(X) 
(i) 

(2) (1) (2) nl 

2: = c\' - a 2 , c 2 = aia 2 - pi 
for i = 3 : k do 

4° = <t x ] - ^t 1 > - ft 

6: for j = 3 : i — 1 do 

c W -c (i_1) -a-c (l_1) -5 2 c (l " 2) 
8: end for 

c i ~ fx i c i-l Pi c i-2 

10: end for 

for i = k + 1 : n do 

» _ „(*-i) _ „. 



12: 



if k > 2 then 

14: 4 = 4 - a iC \ - Pf 

for j = 3 : k do 

(i) (i-l) o2 (»— 2) 

end for 



18: end if 

end for 

20: {Now Cj = cf\ 1 <j <k} 



i 


c (l) 


r Ci) 

c 2 


4° 


1 
2 
3 
4 
5 


(i) 

C\ = -Qi 

(2) (1) 

C\' = Cj — OL2 

(3) (2) 

c l = C l a 3 

(4) (3) 

(5) (4) 

c l = C l «5 


4 2> = maa - /9| 

4 3) = 4 2) -a3C«-A 2 

4 4) = 4 3) -a 4 4 3) -/3 4 2 

4 5) = 4 4) -«s4 4) -/3i 


C W = c <3) _ Q4C <3) _ ^| C W 

4 5) = 4 4) -« 5 4 4) -/3l4 3) 



Fig. 4.1. Coefficients computed by Algorithm 1 for n = 5 and A; = 3. 



4.2. Assumptions for Running Error Bounds. We assume that all matrices 
are real. Error bounds for complex matrices are derived in [261 §6]. In addition, we 
make the following assumptions: 

1. The matrix elements are normalized real floating point numbers. 

. . . ~(i) 

2. The coefficients computed in floating point arithmetic are denoted by c k . 

3. The output from the floating point computation of Algorithms 1 and 2 is 
fl[cfe] = \ In particular, fl[ci] = c\ = — ol\. 

4. The error in the computed coefficients is ei*' so that = and 

= c ( l ] + ef, 2<i<n, 1 < k < n. (4.2) 

5. The operations do not cause underflow or overflow. 

6. The symbol u denotes the unit roundoff, and nu < 1. 

7. Standard error model for real floating point arithmetic [141 §2.2]: 
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If op G {+,— , x,/}, and x and y are real normalized floating point numbers 
so that xo~py does not underflow or overflow, then 



flfxopy] = (xopy)(l + 5) where \S\ <u, (4.3) 

and 

x op y 

H[xopy] = y-j- — where |e| < u. (4-4) 
The following relations are required for the error bounds. 

Lemma 4.1 (Lemma 3.1 and Lemma 3.3 in [14|). Let Si and pi be real numbers, 
1 < i < n, with \5i\ < u and p t — ±1. If nu < 1 then 
1- nr=i(l + ^) W = ! + where 

i/i i nu 

< 7n = ^ • 



a (i + ^)(i + fe ) = i + ^ +fc 

4.3. Running Error Bounds. We derive running error bounds for Algorithm 
1, first for c\, then for £2, and at last for the remaining coefficients Cj, 3 < j < /c. 
The bounds below apply to lines 2, 4, and 14 of Algorithm 1. 
Theorem 4.2 (Error bounds for ). If the assumptions in tU^hold and 



»(*-!) 



2 < i < n, 



|e^| <|e? -1) | + u|c^|, 2<z<n. 



Proo/. The model flU) implies (1 + e^)cf = cf 1} - a, where < u. 

Writing the computed coefficients and c-[ in terms of their errors (|4.2p and 
then simplifying gives 

e f =e H)_ e W a f). 

Hence |e^| < |e^ 1} | + u \c^\. □ 

The bounds below apply to lines 2, 6, and 16 of Algorithm 1. 



Theorem 4.3 (Error bounds for c\ ). If the assumptions in 



hold, 



S (2) 



= fl 



= fl 



fl[aia 2 ] -fl[/? 2 2 ] 



-fl 



CtiC\ 



- fl [ A ? 



3 < « < n, 



then 



S 2 



I "(2) 



c 2 1 < u Mci!2ai| + |/? 2 

4°i < i4 <_1) i + ke£ <_1) i+« + ia 2 i + i4 4) i) +72 k^l- 



s 



Proof. The model (|4.3[) implies for the multiplications 



*(2) - 



a 2 ai(l + 5) -/3f(l + t?) 



where |<5|, | | < u. 



Applying the model (I4.4p to the subtraction gives 



(1 + e)4 2) = a 2 ai(l + <5) - /3 2 2 (1 + 77), |e| < 



Now express c^ 2 " 1 in terms of the errors eif 1 from (I4.2[) and simplify. 

For 3 < i < n, applying model (|4.3|) to the multiplications and the first subtraction 



(2) 



Ji) n 

jives c\ ' = 11 



(0 (i) 

5i - 92 



, where 

gf ee ct'Hl + <5 W ) Oi^Hl + Of), g® = f}f(l + r?W), 

\S^\, \rf 1 } I < u and |0 2 ^ | < 72. Applying the model (|4.4I) to the remaining subtraction 
gives 

(1 + e ('))4 4) = 4^(1 + 5«) - a^-^a + - a 2 (i + „W), 

where |e^| < u. Now express 4* and 4* m terms of their errors (14. 2p to get 

e W = e£-D + _ aie (<-D _ eVoig-V - fir,® ~ e (i) 4°, 

and apply the triangle inequality. □ 

The bounds below apply to lines 8, 10, and 18 of Algorithm 1. 

Theorem 4.4 (Error bounds for Cj \ 3 < j < k). If the assumptions in 
hold, and 



a i c i-l 



M-l) 



«?g(i-2) 
c i-2 



3 < i < k, 



3 < j < k, j + 1 < i < n, 



then 



ef| < M^l + |/3^^ 2) | + u + l#l) +72 I A 

lefl^le^l + la^l + l^^l 



2 4^ 2) l 



Proof. The model (|4.3[) implies for the three multiplications that 



4r 1 1) (i+<5)+/3f4^ 2) (i+^) 



where |<5| < u and |6*2| < 72- Applying model (|4.4[) to the remaining addition gives 

(1 + e)<f = -OitTiHl +&)- ft 2 €"2 2) (l + fc), |e| < u. 
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As in the previous proofs, write c 



s(0 a(<-i) 



and c^i_2^ m terms of their errors 



,(0 



a? «-a) 

Pi e i-2 



2-( 1 -2) 
i L i-2 



Saicll^ -ecf. 



and apply the triangle inequality. 

For k + 1 < i < n, applying (|4.4j) to the two multiplications and the first subtrac- 



fi 



(*) (*) 
Si - 92 



where 



tion gives Cj 

gf ee cf-^il + *«) - ^^(l + $ = pt$:?(i + (f ), 

|<JW| < u and 1^1, \ef\ < 72. Applying model gl} to the remaining subtraction 
gives 

(1 + eW)c« = cf-^il + *«) - ai c^(l + - A 2 C" 2 2) (1 + ft 
where e*- 1 -* < w. Write the computed coefficients in terms of their errors 



and apply the triangle inequality. □ 

We state the bounds when the leading k coefficients of p(X) are computed by 
Algorithm f in floating point arithmetic. 

Corollary 4.5 (Error Bounds for fl(cj), 1 < j < k). If the assumptions ir. 
hold, then 



\ft[cj]-Cj\ < 



l<j<k, 



where H[cj] = c}™ 1 and <pj = \e^\ are given in Theorems \4-£\ [7^1 and \Jm 

5. Nonsymmetric Matrices. In the first stage, La Budde's method [5] reduces 
a real square matrix to upper Hessenberg form H. In the second stage it computes the 
coefficients of the characteristic polynomial of H. We present recursions to compute 
individual coefficients in the second stage of La Budde's method in M 5 . H and derive 
running error bounds in 



5.1. The Algorithm. We present an implementation of the second stage of La 
Budde's method for nonsymmetric matrices. Let 



H = 



(oil hi2 ... 
(3 2 a 2 h 2 3 



V 



Oin ) 



be a real n x n upper Hessenberg matrix with diagonal elements on, subdiagonal 
elements ft, and characteristic polynomial p(X) = det(A7 — H). 

La Budde's method computes the characteristic polynomial of an upper Hes- 
senberg matrix H by successively computing characteristic polynomials of leading 
principal submatrices Hi of order i [5j . Denote the characteristic polynomial of Hi by 
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Pi(X) = det(AI — Hi), 1 < i < n, where p(A) = p n (X). The recursion for computing 
p(X) is [SI (6.57.1)] 

p (A) = 1, pi(A) = A - at 

i-l 

Pi(X) = (A - ai)pi-i(X) - ^ hi-m,i Pi ■ •■Pi-m+l Pi-m-l(A), (5.1) 

m=l 

where 2 < i < n. The recursion for Pi(X) is obtained by developing the determinant 
of XI — Hi along the last row of Hi. Each term in the sum contains an element in the 
last column of Hi and a product of subdiagonal elements. 
As in the symmetric case, we let 

p(X) = A™ + Cl X n - X + ■■■ + Cn-xX + c n 

and 

K (A) = A' + c^A*- 1 + • ■ • + c&A + cf, l<i<n, 

where = c^. Equating like powers of A in (15.11) gives recursions for individual 
coefficients Ck, which are presented as Algorithm 2. In the process, c\,... ,Ck-i are 
also computed. 



Algorithm 2 La Budde's method for upper Hessenberg matrices 
Input: n x n real upper Hessenberg matrix H, index k 
Output: Coefficient Ck of p(X) 



1 

2 
3 
1 
5 
6 
7 
8 
9 
10 
11 
12 
13 
11 
15 
16 
17 
18 



(2) (1) (2) , n 

c\' = c\' — a 2 , c\ = a\a 2 - n 12 p 2 

for i = 3 : k do 

(i) (i-i) 
Ci = q - oti 

for j = 2 : i — 1 do 

(i) (i— 1) (i— 1) v~*7— 2 i r> (i — m — 1) j a a 

C) —c) - aiC}_ a - 2Jm=l A ' • • Pi-m+1 c}_ m _j - /3< • • • A-j+2 

end for 

fi) (i— 1) v-^i— 2 , a a (i—m—1) u a a 

c\ = -a i cl_ 1 - }^ m=1 hi- m ,i Pi--- Pi-m+i c-_ m _i -h u Pi--- p 2 
end for 

for i = k + 1 : n do 

(0 (i-l) 
C\ = C\ -OL{ 

if k > 2 then 

for j = 2 : fc do 

(i) fi— 1) (i— 1) — 2 7 /-) £i (i — m—1) t n 

C) ' = C) - a i C}_ 1 - Em=l fti-m,i Pi--- Pi-m+1 Cj- m -l ~ hi-j+l,i Pi ■ • ■ fii-j+2 

end for 
end if 
end for 

{Now Cj = c[ n \ 1 <j < k} 



For the special case when H is symmetric and tridiagonal, Algorithm 2 reduces 
to Algorithm 1. Figure [57T1 shows an example of the recursions for n = 5 and k = 3. 
Algorithm 2 computes the characteristic polynomial of companion matrices ex- 
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Oil 
„(3) 



f>2 

(4) 



Or, 



4 2) 

c< 3) 



2 
.(4) 
2 

.(5) 
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c (2) 

L 2 

c <3) 

L 2 
.(4) 



aiQ2 

(2) 
«3Ci 

(3) 
(4) 

a 5 q 



/ll2p2 
/l34p4 



„(*) 



c (3) _ 




(2) 






c (4) - 

C 3 ^ 


c (3) 

C 3 


(3) 
— CX4C 2 




- /l24p4p3 


„(5) _ 
C 3 ^ 


c (4) 


(4) 
- Q5C2 


- /l 45 /3 5 4 3) 


- ^35/35/34 



Fig. 5.1. Coefficients Computed by Algorithm 2 when n = 5 and fe = 3. 



actly. To see this, consider the n x n companion matrix of the form 



1 

V 



-c 2 



Algorithm 2 computes c 



for 1 < j < n and 1 < i < n — 1, so that c 



(n) 



Since only trivial arithmetic operations are performed, Algorithm 2 computes the 
characteristic polynomial exactly. 

5.2. Running Error Bounds. We present running error bounds for the coeffi- 
cients of p(X) of a real Hessenberg matrix H. 

The bounds below apply to lines 2, 4, and 11 of Algorithm 2. 

Theorem 5.1 (Error bounds for Cj ). If the assumptions in hold and 



2 < i < n, 



then 



\e?\ < le^l+ul^l 



2 < i < n. 



Proof. The proof is the same as that of Theorem 14.21 □ 

The following bounds apply to line 2 of Algorithm 2, as well as lines 6 and 14 for 
the case j = 2. 

Theorem 5.2 (Error bounds for c 2 )■ If the assumptions in §4^hold, and 



s(2) _ 



fl[aia 2 ]-fl[ftia /3 2 ] 



QijC 



(i-1) 



fl[/l 



3 < i < n, 



then 



,(2) 



| <u (|a 2 ai| + W2| + |4 2) l^ 

+ |a ie ^ 1} | +u (\ct 1} \ + \hi-i,iPi\ + |c 2 l) |) + 72 \ ai cf 



\e?\<\et 1] 
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hold, 



Proof. The proof is the same as that of Theorem l4.3l □ 
The bounds below apply to lines 6,8, and 14 of Algorithm 2. 
THEOREM 5.3 (Error bounds for cj , 3 < j < k). If the assumptions in 



■fl 



fl 



(i-i) 



3 < i < fe, and 
cf =fl II 



a,c 



i-2 



^ ' hi—m,i Pi ' ' ' Pi — m+1 Cj_ m _j + /lli ft "-ft 



r .(»-i) 



•a 



hi-m,i Pi ' ' - ft-m+1 ^i-m-1 ^ ~ ^i—j+l Pi ' ' ' Pi-j+2 



3<j<k, j + l<i<n, then 

i-2 



< l^e^/^ + ^ |/ii-m,i A • • -ft-m+i e^,™-! 1 ' 1 



+7i+l |^i-m,i ft ' ' • ft-m+1 C^y^j 1 

m=2 

+7, (i^-m ft Cfi + N ft ■ • • fti) + « (icfi + Mr^i) 



and 



lefl^le^l + Kef 1 ) 



i-2 



1 I "t~ ^ ' \h-i—m,i ft ' ' ' ft— m+1 e ^-_ m _ 1 | 
m=l 

+7i+l I X! I^-"M ft " ' ft-m+l 4*-m-l I J 
\m=2 / 

+7i (l^i-l.i ft 4^2 2) | + ft ' ••ft-j+2M 

+u (| a f| + |cf 1 )|)+ 72 |a4r i 1 )|. 

Proof The big sum in contains i — 2 products, where each product consists 
of m + 2 numbers. For the m + 1 multiplications in such a product, the model (|4.3j) 
and Lemma |4. II imply 



^m — A hi—m,i ft ' ' ' ft— m+1 ^i — m — 1 — h>i—m,i ft ' ' ' ft— m+1 Cj_ m _i (1 + $m+l), 

where |# m +i| < 7m+i and 1 < m < i — 2. The term hu P% - •• fti is a product of i 
numbers, so that 

= fl [fc M ft • • - A] = ft u ft ••• /3 2 (1 + t -_i), 
13 



where \Oi-i \ < 7i-i- Adding the i — 1 products g m from left to right, so that 

g = &[...&[& \gx + g2]+ga]--'+9i-i], 
gives, again with (|4.3|) . the relation 



i-2 



.9 — /tt-l,t A C^!_ 2 2 (1 + Oi) + ^ hi- m ,i ft ■ ■ ■ ft_m+l 



[i— m— 1) 



q(m) 
7 i+l 



+hu0i~-fa 1 + 



where < 7i+i and < 7;. For the very first term in cj; we get 



a;c. 



(i-i) 



ctiCi—i (1 + <^)j where |<5| < u. Adding this term to <? and using 



model (|4.4I) yields 

-(1 + e)cf } = ^^(l + 5) + ft c£" 2 2) (l + 0;) + h u ft • • ■ ft (1 + 0<) 

i-2 

+ ^2 h*-m,iPi ' ' ' A-m+1 C^^L^ (l + , 
m=2 

where |e| < u. Write the computed coefficients in terms of their errors (|4.2[) 

i-2 

— e i — a i e i-l + 2^ n i~"hk Pi " ' Pi-m+1 e i-m-l + a i c i-l d 
m=\ 

i-2 

I L o 2)/) , i o o ~(i-m— 1) n(m) 

~r"i— 1,1 Pi ^i — 2 i / . ^i— m,i Pi ' ' ' Pi— m+1 Cj_ m _i "i-)-l 

ft • • • ft $i + e a 
and then apply the triangle inequality 

(i) 

For j + 1 < i, contains the ; 
first subtraction. Model (|4.3I) implies 



m=2 



~(i) . . . ~(i — 1) 

For j + 1 < i, Cj- contains the additional term , which is involved in the 



fl 



4-)-fi 



_(i-i) 



= <7 



- 1} (l + ^-^ti 1 ' (l + ^) 



where \ < u and I < 72- From this we subtract g which is computed as in the 
case j = i. □ 

Finally we can state bounds when the leading k coefficients of p{\) are computed 
by Algorithm 2 in floating point arithmetic. 

Corollary 5.4 (Error Bounds for fi(cj), 1 < j < k). If the assumptions in §^.g| 
hold, then 



|n[ Cj ]-c,| < Pj , l<j<k, 



where fl[cj] = Cj and pj = |e^™' l | are given in Theorems \5.1[ \5.S\ and \5.SX 

Potential instability of La Budde 's method. The running error bounds reflect the 
potential instability of La Budde's method. The coefficient Cj is computed from the 
preceding coefficients c^ 1 ^ c^~^ +1 ^ . La Budde's method can produce inaccurate 
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results for c) , if the magnitudes of preceding coefficients are very large compared to 

(i) ... . (i) 

Cj so that catastrophic cancellation occurs in the computation of Cj . This means 
the error in the computed coefficient Cj can be large if the preceding coefficients in 
the characteristic polynomials of the leading principal submatrices are larger than Cj . 

It may be that the instability of La Budde's method is related to the illcondition- 
ing of the coefficients. Unfortunately we were not able to show this connection. 

6. Overall Error Bounds. We present first order error bounds for both stages 
of La Budde's method. The bounds take into the account the error from the reduction 
to Hessenberg (or tridiagonal) form in the first stage, as well as the roundoff error from 
the computation of the characteristic polynomial of the Hessenberg (or tridiagonal) 
matrix in the second stage. We derive bounds for symmetric matrices in Q6.l[ and for 
nonsymmetric matrices in H6.2\ 

6.1. Symmetric Matrices. This bound combines the errors from the reduction 
of a symmetric matrix A to tridiagonal form T with the roundoff error from Algorithm 
1. 

Let T = T + E be the tridiagonal matrix computed in floating point arithmetic 
by applying Householder similarity transformations to the symmetric matrix A. From 
[TU1 §8.3.1.] follows that for some small constant v\ > one can bound the error in 
the Frobenius norm by 

\\E\\ F < Vl n 2 \\A\\ F u. (6.1) 

The backward error E can be viewed as a matrix perturbation. This means we 
need to incorporate the sensitivity of the coefficients Cj to changes E in the matrix. 
The condition numbers that quantify this sensitivity can be expressed in terms of 
elementary symmetric functions of the singular values |19j . Let a\ > . . . > a n be the 
singular values of A, and denote by 

s = 1, Sj = ^2 <Tii"-<Ti j > 1 < i < n, 

l<ii<---<ij<n 

the jth elementary symmetric function in all n singular values. 

Theorem 6.1 (Symmetric Matrices). If the assumptions in §^.l?| hold, A is real 
symmetric with \\A\\p < \/{vin 2 u) for the constant v\ in \6. Cj are the coefficients 
of the characteristic polynomial of T , then 

I fl[c,-] - Cj| < (n - j + 1) Sj-i ^in 2 ||A|| F u + <j)j+0 (u 2 ) , 1 < j < k, 

where <fij are the running error bounds from Corollary \4-5\ 
Proof. The triangle inequality implies 

I fl[cj] - Cj \ < | fi[cj] -£j\ + \cj -Cj\. 

Applying Corollary 15.41 to the first term gives | fl[cj] — Cj \ < (f>j. 

Now we bound the second term \cj — Cj\, and use the fact that A and T have the 
same singular values. If \\E\\2 < 1 then the absolute first order perturbation bound 
[191 Remark 3.6] applied to T and T + E gives 

\cj-Cj\ < (n-j + ^Sj-x \\E\\ 2 u + 0(\\E\\ 2 2 ) , 1 <j< k. 
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From (jlTTj) follows [[_B|| 2 < \\E\\ F < Vi<n 2 \\A\\ F u. Hence we need \\A\\ F < \j{v x n 2 u) 
to apply the above perturbation bound. □ 

Theorem 16.11 suggests two sources for the error in the computed coefficients fl|cj]: 
the sensitivity of cj to perturbations in the matrix, and the roundoff error pj in- 
troduced by Algorithm 1. The sensitivity of Cj to perturbations in the matrix is 
represented by the first order condition number (n — j + l)sj_i, which amplifies the 
error j/in 2 ||A||i? u from the reduction to tridiagonal form. 

6.2. Nonsymmetric Matrices. This bound combines the errors from the re- 
duction of a nonsymmetric matrix A to upper Hessenberg form H with the roundoff 
error from Algorithm 2. 

Let H = H + E be the upper Hessenberg matrix computed in floating point 
arithmetic by applying Householder similarity transformations to A. From [10l §7.4.3] 
follows that for some small constant v 2 > 

\\E\\ F < ^n 2 \\A\\ F u. (6.2) 

The polynomial coefficients of nonsymmetric matrices are more sensitive to changes 
in the matrix than those of symmetric matrices. The sensitivity is a function of only 
the largest singular values, rather than all singular values [15]. We define 

4 = 1, s/2i = On • • • < j<ri ■ ■ ■ <Tj-i, 1 < j < n, 

l<h<---<ij-i<j 

which is the (j — l)st elementary symmetric function in only the j largest singular 
values. 

Theorem 6.2 (Nonsymmetric Matrices) . If the assumptions in §4-2\ hold, \\A\\ F < 
l/(u27i 2 u) for the constant V2 in and Cj are the coefficients of the characteristic 

polynomial of H , then 

Iflfe] - Cj\ < Q 4-1 v ^ 2 \\ a \\f u + Pj + O (u 2 ) , l<j<k, 

where pj are the running error bounds from Corollary \5.J\ 

Proof. The proof is similar to that of Theorem l6.ll The triangle inequality implies 

&[Cj] -Cj\<\ &[Cj] -Cj\ + \Cj -Cj\. 

Applying Corollary [53] to the first term gives | fl.[£j] — Cj\ < Pj- 

Now we bound the second term \5j — Cj\, and use the fact that A and H have the 
same singular values. If ||-E||2 < 1 then the absolute first order perturbation bound 
[151 Remark 3.4] applied to H and H + E gives 




1 < j < k. 



From (JO) follows |[^[[ 2 < \\E\\ F < v 2 n 2 \\A\\ F u. Hence we need \\A\\ F < l/{v 2 n 2 u) 
to apply the above perturbation bound. □ 

As in the symmetric case, there are two sources for the error in the computed 
coefficients fl[cj]: the sensitivity of Cj to perturbations in the matrix, and the roundoff 
error pj introduced by Algorithm 2. The sensitivity of Cj to perturbations in the 
matrix is represented by the first order condition number (ps^^, which amplifies 
the error ^2^ 2 ||^4||f u from the reduction to Hessenberg form. 
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Fig. 7.1. Forsythe Matrix. Lower (blue) curve: Absolute errors \c£ b2 (F) — ci c (F)\ of the 
coefficients computed by Algorithm 2. Upper (red) curve: Running error bounds pk from Corollary 




20 40 60 80 100 120 140 160 180 200 



Fig. 7.2. Forsythe Matrix. Coefficients c^ oly (F) computed by poly. The exact coefficients 
c k = 0, 1 < k < 199. 



7. Numerical Experiments. We compare the accuracy of Algorithms 1 and 2 
to MATLAB's poly function, and demonstrate the performance of the running error 
bounds from Corollaries 14.51 and 15.41 The experiments illustrate that Algorithms 1 
and 2 tend to be more accurate than poly, and sometimes substantially so, especially 
when the matrices are indefinite or nonsymmetric. 

We do not present plots for the overall error bounds in Theorems 16.11 and 16.21 
because they turned out to be much more pessimistic than expected. We conjecture 
that the errors from the reduction to Hessenberg form have a particular structure that 
is not captured by the condition numbers. 



The coefficients computed with Algorithms 1 and 2 are denoted by 



algl 



algl 



respectively, while the coefficients computed by poly are denoted by 



and 

„poly 



Furthermore, we distinguish the characteristic polynomials of different matrices by 
using Cfc(A) for the fcth coefficient of the characteristic polynomial of the matrix X . 
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Fig. 7.3. Hansen's Matrix. Upper (blue) curve: Relative errors \c^° ly (H) — Ck(H)\/\ck(H)\ 

of coefficients computed by poly. Lower (red) curve: Relative errors \c£ 9l (H) — Ck(H)\/\ck(H)\ of 
coefficients computed by Algorithm 1. 



7.1. The Forsythe Matrix. This example illustrates that Algorithm 2 can 
compute the coefficients of a highly nonsymmetric matrix more accurately than poly, 
and that the running error bounds from Corollary 15 .41 approximate the roundoff error 
from Algorithm 2 well. 

We choose a n x n Forsythe matrix, which is a perturbed Jordan block of the 
form 



/0 1 



\ 

o 1 
o / 



where v = 10 



-10 



(7.1) 



with characteristic polynomial p(X) = A" — v. Then we perform an orthogonal simi- 
larity transformation F\ = QF2Q T , where Q is an orthogonal matrix obtained from 
the QR decomposition of a random matrix. The orthogonal similarity transformation 
to upper Hessenberg form F is produced by F = hess(Fi). 

We applied Algorithm 2 to a matrix F of order n = 200. Figure 17.11 shows that 
Algorithm 2 produces absolute errors of about 10~ 15 , and that the running error 
bounds from Corollary 15.41 approximate the roundoff error from Algorithm 2 well. In 
contrast, the absolute errors produced by poly are huge, as Figure [7T21 shows. 

7.2. Hansen's Matrix. This example illustrates that Algorithm 1 can com- 
pute the characteristic polynomial of a symmetric positive definite matrix to machine 
precision. 

Hansen's matrix |12[ p 107] is a rank one perturbation of a n X n symmetric 
tridiagonal Toeplitz matrix, 



H 



f 1 1 

-1 2 

V 



-1 

2/ 
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Fig. 7.4. Hansen's Matrix. Lower (blue) curve: 
efficients computed by Algorithm 1. Upper (red) curve. 



Absolute errors \c^( gl (H) — Cfc| in the co- 
Running error bounds <j>k from Corollary 



"I ±±± + + . +++.+ h ++ ++, +++, + ++ +++++++ 

10 20 30 40 50 60 70 80 90 100 
even k 



Fig. 7.5. Symmetric Indefinite Tridiagonal Toeplitz Matrix. Upper (blue) curve: Relative 
errors |c^° V (T) — c^ (T)\/\ck(T)\ in the coefficients computed by poly for even k. Lower (red) curve: 
Relative errors |c^ 9l (T) — cj.(T)|/|cfc(T)| in the coefficients computed by Algorithm 1 for even k. 

Hansen's matrix is positive definite, and the coefficients of its characteristic polyno- 
mial are 

c n . k+1 (H) = ^ + * " ^ , l < k < n . 

Figure [731 illustrates for n — 200 that the Algorithm 1 computes the coefficients to 
machine precision, and that later coefficients have higher relative accuracy than those 
computed by poly. With regard to absolute errors, Figure 17.41 indicates that the 
running error bounds (f>j from Corollary 14.51 reflect the trend of the errors, but the 
bounds become more and more pessimistic for larger k. 

7.3. Symmetric Indefinite Toeplitz Matrix. This example illustrates that 
Algorithm 1 can compute the characteristic polynomial of a symmetric indefinite 
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Fig. 7.6. Symmetric Indefinite Tridiagonal Toeplitz Matrix. Lower (blue) curve: Absolute 
errors |c^ !sl (T) — Cfe(T)| in the coefficients computed by Algorithm 1. Middle (red) curve: Running 
error bounds <j>j from Corollary \4-5\ Upper (green) curve: Absolute errors \c^ oly (T) — Cfc(T)| in the 
coefficients computed by poly. 



matrix to high relative accuracy, and that the running error bounds in Corollary 14.51 
capture the absolute error well. 

The matrix isanxn symmetric indefinite tridiagonal Toeplitz matrix 

/0 100 \ 

/ 100 '•• '■■ 

100 
V 100 / 

where the coefficients with index are zero, i.e. C2j— i(T) = for j > 1. 

For n = 100 we obtained the exact coefficients Cfc(T) with sym2poly (poly (sym(T) ) ) 
from MATLAB's symbolic toolbox. Algorithm 1 computes the coefficients with odd 
index exactly, i.e. c 2 j_i(T) = for j > 1 and 1 < i < n. In contrast, as Figure 17761 
shows, the coefficients computed by poly can have magnitudes as large 10 185 . 

Figure [7751 illustrates that Algorithm 1 computes the coefficients C2j(T) with even 
index to machine precision, while the coefficients computed with poly have relative 
errors that are many magnitudes larger. Figure [7761 also shows that the running error 
bounds approximate the true absolute error very well. What is not visible in Figure 
17.61 but what one can show from Theorems 14.21 and 14.31 is that </>2j-i = 0. Hence the 
running error bounds recognize that C2j-\ are computed exactly. 

7.4. Frank Matrix. This example shows that Algorithm 2 is at least as ac- 
curate, if not more accurate than poly for matrices with ill conditioned polynomial 
coefficients. 

The Frank matrix U is an upper Hessenberg matrix with determinant 1 from 
MATLAB's gallery command of test matrices. The coefficients of the characteristic 
polynomial appear in pairs, in the sense that Ck(U) — c n -k{U). For a Frank matrix of 
order n = 50, we used MATLAB's toolbox to determine the exact coefficients Ck(U) 
with the command sym2poly (poly (sym(U) ) ) . Figure 17771 illustrates that Algorithm 
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Fig. 7.7. Prank Matrix. Upper (blue) curve: Absolute errors \c^° ly (U) — Ck(U)\ in the coeffi- 
cients computed by poly. Lower (red) curve: Absolute errors \c£ a2 (U) — c^{U)\ in the coefficients 
computed by Algorithm 2. 




10 12 14 16 18 20 22 24 26 



Fig. 7.8. Frank Matrix. Relative errors |c? !ff2 (C/) — Ck(U)\/\ck(U )| in the first 25 coefficients 
computed by Algorithm 2. 



2 computes the coefficients at least as accurately as poly. In fact, as seen in Figure 
I7.8[ Algorithm 2 computes the first 20 coefficients to high relative accuracy. 



7.5. Chow Matrix. This example illustrates that the errors in the reduction to 
Hcssenberg form can be amplified substantially when the coefficients of the charac- 
teristic polynomial are illconditioned. 

The Chow matrix is a matrix that is Toeplitz as well as lower Hessenberg from 
MATLAB's gallery command of test matrices. Our version of the transposed Chow 
matrix is an upper Hessenberg matrix with powers of 2 in the leading row and trailing 
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Fig. 7.9. Transposed Chow Matrix. Upper (blue) curve: Relative errors \c^° ly (C T ) — 
c fc(C T )l/l c fc(C T )l * n the coefficients computed by poly. Lower (red) curve: Relative errors 
\c£ a2 {C T ) — Ck(C T )\/\ck(C T )\ in the coefficients computed by Algorithm 2. 



5 10 15 20 25 30 35 40 45 50 



Fig. 7.10. Chow Matrix. Blue curve: Relative errors \c%° y (C) — c k (C)\/\ck-(C)\ in the coeffi- 
cients computed by poly. Red curve: Relative errors |c^ is2 (C) — Cfc(C)|/|cj.(C)| in the coefficients 
computed by Algorithm 2. The two curves are virtually indistinguishable. 



column, 

/3 4 ... 2™\ 

C r 1 '• '• : . 

'•• 3 4 
V 1 3/ 

As before, we computed the exact coefficients with MATLAB's symbolic toolbox. 
Figure lTl?! illustrates that Algorithm 2 computes all coefficients Ck(C T ) to high relative 
accuracy for n = 50. In contrast, the relative accuracy of the coefficients computed 
by poly deteriorates markedly as k becomes larger. 

However, if we compute instead the characteristic polynomial of C, then a pre- 
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liminary reduction to upper Hessenberg form is necessary. Figure [7. 101 illustrates that 
the computed coefficients have hardly any relative accuracy to speak of, and only the 
trailing coefficients have about 1 significant digit. The loss of accuracy occurs because 
the errors in the reduction to Hessenberg form are amplified by the condition numbers 
of the coefficients, as the absolute bound in Theorem 16.21 suggests. Unfortunately, in 
this case, the condition numbers in Theorem 16.21 are too pessimistic to predict the 
absolute error of Algorithm 2. 

Acknowledgements. We thank Dean Lee for helpful discussions. 
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